#' Export a Formattable as PNG, PDF, or JPEG
#'
#' @param f A formattable.
#' @param file Export path with extension .png, .pdf, or .jpeg.
#' @param width Width specification of the html widget being exported.
#' @param height Height specification of the html widget being exported.
#' @param background Background color specification.
#' @param delay Time to wait before taking webshot, in seconds.
#'
#' @importFrom formattable as.htmlwidget
#' @importFrom htmltools html_print
#' @importFrom webshot webshot
#'
#' @export
export_formattable <- function(f, file, width = "100%", height = NULL, 
                               background = "white", delay = 0.2)
{
  w <- as.htmlwidget(f, width = width, height = height)
  path <- html_print(w, background = background, viewer = NULL)
  url <- paste0("file:///", gsub("\\\\", "/", normalizePath(path)))
  webshot(url,
          file = file,
          selector = ".formattable_widget",
          delay = delay)
}

###############################################################################
################### MULTIPLE HISTOGRAMS IN ONE FIGURE FUNCTION ################
###############################################################################

## Functions you'll need

plot_multi_histogram <- function(df, feature, label_column) {
  plt <- ggplot(df, aes(x = eval(parse(text=feature)), fill = eval(parse(text = label_column)))) +
    geom_histogram(alpha = 0.7, position="identity", 
                   aes(y = ..density..), color = "black") +
    geom_density(alpha = 0.7) +
    geom_vline(aes(xintercept = 0), color = "black", linetype = "dashed", size = 1) +
    labs(x = feature, y = "Density") + 
    theme_custom()    
  plt + guides(fill = guide_legend(title = label_column)) +
    ggtitle("Effect Size Distribution by Tax Type")
}


###############################################################################
################## FUNCTION TO TAKE RESULTS FROM exBounds #####################
###############################################################################

## Input needs to have the following form:
## List of length 2, first object a data frame of model statistics
## second object coefficients and associated se, t, and p

exStats <- function(resList = NULL){ # Needs to have particular form: 
  
  
  modStats <- lapply(resList, FUN = function(x){x[["mod_stats"]]})
  dfMod <- as.data.frame(round(do.call("rbind", modStats), digits = 4))
  dfEst <- do.call("rbind.fill", 
                   lapply(resList, FUN = function(x){x[["mod_coefs"]]}))

  ## Format the data frame
  
  colnames(dfMod) <- c("diff_coef", "diff_se", "diff_t", "diff_p",
                       "rsq", "adj_rsq", "D_stat", "nobs", "prop_nonmiss", 
                       "num_countries", "t_min", "t_max", "t_bar") # make sure these names match up with order in newStats
  dfMod[1:length(colnames(dfMod))] <- round(sapply(dfMod[1:length(colnames(dfMod))], as.numeric), digits = 4)
  
  ## Unpack the data frame dfEst
  
  listEst <- list()
  qoi <- c("beta", "se", "t", "p")
  for(i in 1:length(qoi)){
    listEst[[qoi[i]]] <- dfEst[seq(from = i, to = nrow(dfEst), by = length(qoi)), ]
  }
  
  
  ## Average and median estimates 
  
  (dfCoefsAvg <- round(simplify2array(lapply(listEst, FUN = function(x){apply(x, 2, mean, na.rm = T)})), digits = 4))
  
  ## Some summary information about the coefficients from the models
  
  propNeg <- apply(listEst[[1]], 2, FUN = function(x){
    mean(ifelse(x < 0, 1, 0), na.rm = T)})
  propSig <- apply(listEst[[4]], 2, FUN = function(x){
    mean(ifelse(x < 0.1, 1, 0), na.rm = T)})
  posCoefs <- apply(listEst[[1]], 2, 
                    FUN = function(x){ifelse(x > 0, 1, -1)})
  sigCoefs <- apply(listEst[[4]], 2, 
                    FUN = function(x){ifelse(x < 0.1, 1, 0)})
  sigNeg <- apply(posCoefs*sigCoefs, 2, 
                  FUN = function(x){mean(ifelse(x == -1, 1, 0), na.rm = T)})
  sigPos <- apply(posCoefs*sigCoefs, 2, 
                  FUN = function(x){mean(ifelse(x == 1, 1, 0), na.rm = T)})
  dfSigAvg <- round(data.frame(prop_neg = propNeg, 
                               prop_pos = 1 - propNeg, 
                               prop_sigpos = sigPos,
                               prop_signeg = sigNeg), digits = 4)
  
  ## Weighted average
  
  #weights <- dfMod$adj_rsq/sum(dfMod$adj_rsq, na.rm = T) # Alternative weighting that only weights by r-squared
  weights <- c(dfMod$prop_nonmiss*dfMod$adj_rsq/sum(dfMod$prop_nonmiss*dfMod$adj_rsq)) # Weights by product of r squared and % non-missing
  
  (dfCoefsWeighted <- round(simplify2array(lapply(listEst, FUN = function(x){
    apply(x, 2, FUN = function(x){weighted.mean(x, w = weights, na.rm = T)})})), digits = 4))
  
  
  ## Some summary information about the coefficients from the models
  
  propNeg <- apply(listEst[[1]], 2, FUN = function(x){
    weighted.mean(ifelse(x < 0, 1, 0), w = weights, na.rm = T)})
  propSig <- apply(listEst[[4]], 2, FUN = function(x){
    weighted.mean(ifelse(x < 0.1, 1, 0), w = weights, na.rm = T)})
  posCoefs <- apply(listEst[[1]], 2, 
                    FUN = function(x){ifelse(x > 0, 1, -1)})
  sigCoefs <- apply(listEst[[4]], 2, 
                    FUN = function(x){ifelse(x < 0.1, 1, 0)})
  sigNeg <- apply(posCoefs*sigCoefs, 2, 
                  FUN = function(x){weighted.mean(ifelse(x == -1, 1, 0), 
                                                  w = weights, na.rm = T)})
  sigPos <- apply(posCoefs*sigCoefs, 2, 
                  FUN = function(x){weighted.mean(ifelse(x == 1, 1, 0), 
                                                  w = weights, na.rm = T)})
  dfSigWeighted <- round(data.frame(prop_neg = propNeg, 
                                    prop_pos = 1 - propNeg, 
                                    prop_sigpos = sigPos,
                                    prop_signeg = sigNeg), digits = 4)
  
  ## Add this to higher level list
  
  dfList <- list(output = listEst, 
                 model = dfMod, # This gives model summary statistics: rsq, adjrsq, artest, number of observations, proportion of non-missing observations, number of countries in sample, minimum number of time periods per country, maximum, average
                 sig_avg = dfSigAvg,
                 sig_weighted = dfSigWeighted, 
                 coefs_avg = dfCoefsAvg,
                 coefs_weighted = dfCoefsWeighted)
  
  return(dfList)
}

###############################################################################
################# WRAPPER THAT WILL DO ONE-RUN OF EXTREME BOUNDS ##############
###############################################################################

exBounds <- function(data = NULL, # Has to be p.series object that is balanced
                    formula = NULL, 
                    conductAR = T, # If TRUE, conducts modified Bhargava Panel DW-test
                    minYears = 2, # Minimum number of years that a country can be in the dataset after accounting for missingness
                    seType = "cluster", # Can be "classic", "cluster", or "robust" 
                    levCluster = "time", # Can be "time", or "unit"
                    whichVar = "within", # Can be "within", "random", "ht", "between", "pooling", "fd"
                    whichFE = "individual"){ # Can be "individual", "time", "twoways" 
  
  dfUse <- data
  useEq <- as.formula(formula)
  useVars <- all.vars(useEq)
  panelInfo <- tapply(dfUse$year, dfUse$ccode, length)
  rmObs <- which(dfUse$ccode %in% as.numeric(names(panelInfo[panelInfo < minYears]))) # Remove any countries that only show up once, can't take a mean there
  if(length(rmObs > 0)){dfUse <- dfUse[-rmObs, ]} # Remove any rows that don't fit criteria above
  useInfo <- tapply(dfUse$year, dfUse$ccode, length)
  
  ## Estimate model
  mod <- plm(useEq, model = whichVar, effect = whichFE, data = dfUse, 
             index = c("ccode", "year"))
  
  ## Format output and create vcov for coefficient test
  if(conductAR){arStat <- pbnftest(mod)$statistic}else{arStat <- NA} # Note that because panel is unbalanced test is modified Durbin-Watson test, recommendation is statistic < 2 suggests no positive serial correlation
  if(seType == "classic"){
    vcovUse <- vcov(mod)
    estMat <- as.data.frame(t(summary(mod)$coefficients[, 1:4]))}
  if(seType == "cluster"){
    vcovUse <- vcovHC(mod, cluster = levCluster)
    estMat <- as.data.frame(t(coeftest(mod, vcov. = vcovUse)[, 1:4]))}
  if(seType == "robust"){
    vcovUse <- vcovHC(mod)
    estMat <- as.data.frame(t(coeftest(mod, vcov. = vcovUse)[, 1:4]))}
  
  rownames(estMat) <- c("beta", "se", "t", "p")
  
  ## Compare coefficients on direct and indirect taxation
  
  coefs <- summary(mod)$coefficients
  betaDiff <- coefs[rownames(coefs) == directVar, 1] - coefs[rownames(coefs) == indirectVar, 1]
  locDirect <- which(rownames(vcovUse) == directVar)
  locIndirect <- which(rownames(vcovUse) == indirectVar)
  seDiff <- sqrt(vcovUse[locDirect, locDirect] 
                 + vcovUse[locIndirect, locIndirect]
                 - 2*vcovUse[locDirect, locIndirect])
  newStats <- c(betaDiff, # difference between indirect and direct taxation coefficients
                seDiff, # se of difference
                betaDiff/seDiff, # t of difference
                2*(1 - pnorm(abs(betaDiff/seDiff))), # p of difference (note this is two-sided)
                summary(mod)$r.squared, # r-squared
                arStat, # modified Durbin-Watson test
                nrow(model.frame(mod)),
                nrow(model.frame(mod))/nrow(dfTax), # This is the % of non-missing data
                length(unique(dfUse$ccode)), 
                min(useInfo, na.rm = T), max(useInfo, na.rm = T), 
                mean(useInfo, na.rm = T))
  names(newStats) <- c("diff_coef", "diff_se", "diff_t", "diff_p",
                       "rsq", "adj_rsq", "D_stat", "nobs", "prop_nonmiss", 
                       "num_countries", "t_min", "t_max", "t_bar") # make sure these names match up with order in newStats
  return(list(mod_stats = newStats, 
              mod_coefs = estMat, 
              mod_formula = formUse))
  
} 

###################################################################
################### CREATE LATEX TABLE FROM DATAFRAME #############
###################################################################

## Shape should be coefficients in the rows, QOIs and model
## differentiating variable in the columns

prepTexreg <- function(data = NULL, 
                       names.col = "term", # Should be names of coefficients
                       coef.col = "estimate", 
                       p.col = "p.value", 
                       se.col = "std.error",
                       lci.col = "conf.low", 
                       uci.col = "conf.high",
                       model.col = "outcome", 
                       alpha = 0.05){ # Names of goodness of fit statistics
                       
  require(multiwayvcov)
  df <- data
  texreg_obj <- list() # This is where you store the texreg objects
  arg_names <- c(names.col, coef.col, p.col, se.col, 
                 lci.col, uci.col, model.col)
  arg_na <- arg_names[!arg_names %in% ls(df)]
  
  if(any(c(lci.col, uci.col) %in% arg_na)){
    df[, lci.col] <- df[, coef.col] - qnorm(1 - alpha/2)*df[, se.col]
    df[, uci.col] <- df[, coef.col] + qnorm(1 - alpha/2)*df[, se.col]
    warning("No confidence intervals supplied or default column 
    names do not match those in data frame. Creating symmetric confidence intervals from
            standard error.")}
  
  if(class(df[, model.col]) == "character")
    mod_names <- unique(df[, model.col])
  if(class(df[, model.col]) == "factor")
    mod_names <- levels(df[, model.col])
  
  for(i in 1:length(mod_names)){
    df_sub <- df[df[, model.col] == mod_names[i], ]
    
    if(any(c(lci.col, uci.col) %in% arg_na))
      texreg_obj[[i]] <- createTexreg(coef = df_sub[, coef.col], 
                                      pvalues = df_sub[, p.col], 
                                      coef.names = df_sub[, names.col], 
                                      se = df_sub[, se.col], 
                                      model.name = mod_names[i])
    
    if(!all(c(lci.col, uci.col) %in% arg_na))
      texreg_obj[[i]] <- createTexreg(coef = df_sub[, coef.col], 
                                      pvalues = df_sub[, p.col], 
                                      coef.names = df_sub[, names.col], 
                                      se = df_sub[, se.col], 
                                      ci.low = df_sub[, lci.col], 
                                      ci.up = df_sub[, uci.col], 
                                      model.name = mod_names[i])
    
      
    }
  
  return(list(obj = texreg_obj, 
              mod_names = mod_names))
  
}
  
###################################################################
################### CUSTOM GGPLOT FUNCTION ########################
###################################################################

theme_custom <- function(legend_justification = c(1, 0), 
                         legend_position = c(0.95, 0.8), 
                         legend_box = "vertical", 
                         title_size = 12,
                         title_position = 0.5, # 0 = left, 1 = right
                         subtitle_position = 0.5,
                         axis_text_size = 10, 
                         axis_title_size = 12,
                         axis_border_size = 0,
                         text_size = 12,
                         panel_border_size = 0.7){
  theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          panel.background = element_blank(),
          title = element_text(size = title_size),
          axis.text = element_text(size = axis_text_size),
          axis.title = element_text(size = axis_title_size),
          axis.line = element_line(size = axis_border_size),
          legend.justification = legend_justification, 
          legend.position = legend_position,
          legend.box = legend_box,
          panel.border = element_rect(colour = "black", fill = NA, size = panel_border_size), 
          panel.spacing = unit(0, "mm"), 
          plot.title = element_text(hjust = title_position), 
          plot.subtitle = element_text(hjust = subtitle_position))  
}

##################################################################################################
############### MAKE EQUATION AS CHARACTER STRING FROM DVs and IVs ###############################
##################################################################################################

make_equation <- function(dvs, ivs, interactions = c()){
  return(paste(dvs, "~", paste(c(ivs, interactions), sep = "", collapse = " + ")))
}


###################################################################
################### QUICK REGRESSION + VISUALIZATION ##############
###################################################################

quick_lm <- function(dep.vars = NULL, 
                     ind.vars = NULL, # note that the order these variables are in will control the order they appear in the plot
                     data = NULL, 
                     alpha = 0.05,
                     dep.names = NULL,
                     data.names = NULL,
                     stand.vars = NULL, # Set this to "non-binary" if you want all non-binary, non-factor/character class variables standardized
                     stand.mult = 2, # This is the multiplier for the standard deviation scaling on plots. Default is 2 to make it more comparable to binary predictors
                     inc.mod = F,
                     inc.plot = F,
                     print.plot = F,
                     keep.vars.tab = NULL, 
                     keep.dep.tab = NULL,
                     keep.data.tab = NULL,
                     keep.vars.plot = NULL, 
                     keep.dep.plot = NULL, 
                     keep.data.plot = NULL,
                     title.plot = c(""), 
                     x.size = 14, 
                     y.size = 12, 
                     strip.size = 14,
                     subset.by = "dataset", # Can be either "dataset" or "outcome"
                     ...){
  
  ################ PACKAGES AND CHECKS #########################
  
  packages.needed <- c("lmtest", "clubSandwich", 
                       "dotwhisker", "broom", "estimatr", 
                       "dplyr")
  if(!all(packages.needed %in% rownames(installed.packages()))){
    stop(paste("Need to install the following packages:",
               packages.needed[!packages.needed %in% rownames(installed.packages())]))}
  lapply(packages.needed, require, character.only = TRUE)
  
  if(class(data) != "list"){data <- list(data)}
  if(any(dep.vars %in% ind.vars)){stop("At least one dependent variable is also be included as a predictor")}
  if(is.null(dep.names)){dep.names <- dep.vars}
  if(is.null(data.names)){data.names <- paste("dataset", 1:length(data), sep = "")}
  if(is.null(subset.by) & length(dep.vars) == 1 & length(data) > 1){
    subset.by <- "outcome"}
  if(is.null(subset.by) & length(dep.vars) > 1 & length(data) == 1){
    subset.by <- "dataset"}
  if(is.null(keep.vars.tab)){keep.vars.tab <- c(ind.vars, "(Intercept)")}
  if(is.null(keep.dep.tab)){keep.dep.tab <- dep.vars}
  if(is.null(keep.data.tab)){keep.data.tab <- data.names}
  if(is.null(keep.vars.plot)){keep.vars.plot <- c(ind.vars, "(Intercept)")}
  if(is.null(keep.dep.plot)){keep.dep.plot <- dep.vars}
  if(is.null(keep.data.plot)){keep.data.plot <- data.names}
  
  ################ ESTIMATION ########################
  
  
  est.eq <- make_equation(dep.vars, ind.vars)
  data.tab <- info.tab <- c()
  get.stats <- c("N", "k", "r.squared", "adj.r.squared")
  mod.data <- list() # This will store model objects
  
  ## Loop over datasets
  
  for(i in 1:length(data)){
  
    ## Check to make sure variable is in each dataset   
    if(!all(c(ind.vars, dep.vars) %in% colnames(data[[i]]))){
      all.vars <- c(ind.vars, dep.vars)
      miss.vars <- all.vars[!all.vars %in% colnames(data[[i]])]
      if(length(miss.vars > 0)){
        stop(paste("The following variables are not present in", data.names[i], ":", 
                   paste(miss.vars, collapse = ", ")))
      }
    }
    
    ## Standardize according to stand.vars argument
    if(all(stand.vars == "non-binary") & !is.null(stand.vars)){
      all.vars <- c(ind.vars, dep.vars)
      change.vars <- all.vars[unlist(lapply(data[[i]][, all.vars], 
             FUN = function(x){
               ifelse(length(unique(x[!is.na(x)])) > 2 & class(x) %in% c("numeric", "integer"), T, F)}))]
      data[[i]][, change.vars] <- scale(data[[i]][, change.vars], 
                                        scale = stand.mult*apply(data[[i]][, change.vars], 2, sd, na.rm = T))
    }
    
    if(all(stand.vars != "non-binary") & !is.null(stand.vars)){
      all.vars <- c(ind.vars, dep.vars)
      add.vars <- stand.vars[!stand.vars %in% c(ind.vars, dep.vars)]
      if(length(add.vars > 0)){
        stop(paste("The following variables are set for standardization but are not in the dependent or independent variable vector", 
                   add.vars))}
      
      change.vars <- stand.vars[unlist(lapply(data[[i]][, stand.vars], 
                                           FUN = function(x){
                                             ifelse(length(unique(x[!is.na(x)])) > 2 & class(x) %in% c("numeric", "integer"), T, F)}))]
      data[[i]][, change.vars] <- scale(data[[i]][, change.vars], 
                                        scale = stand.mult*apply(data[[i]][, change.vars], 2, sd, na.rm = T))
    }
    
    ## Build necessary objects
    coef.tab <- stats.tab <- c()
    mod.dep <- list()
    ## Loop over dependent variables
    for(j in 1:length(dep.vars)){
      main.args <- list(formula = as.formula(est.eq[j]), 
                        data = data[[i]], 
                        alpha = alpha)
      add.args <- list(...)
      mod <- mod.dep[[j]] <- do.call(lm_robust, args = c(main.args, add.args))
      tmp.tab <- tidy(mod)
      stats.tab <- rbind(stats.tab, c(data.names[i], dep.names[j], 
                                      mod[get.stats], mod$fstatistic[1]))
      coef.tab <- rbind(coef.tab, tmp.tab)                  
      
    }
    names(mod.dep) <- dep.vars
    mod.data[[i]] <- mod.dep
    
    ## Clean up the info and tab objects a bit and stack them
    stats.tab <- as.data.frame(stats.tab)
    colnames(stats.tab) <- c("dataset", "outcome", get.stats, "fstatistic")
    info.tab <- rbind(info.tab, stats.tab)
    coef.tab$dataset <- data.names[i]
    data.tab <- rbind(data.tab, coef.tab)
  }
  names(mod.data) <- data.names
  
  colnames(data.tab) <- c("term", "estimate", "std.error", "statistic", 
                          "p.value", "conf.low", "conf.high", "df", 
                          "outcome", "dataset")
  data.tab <- data.tab[,  c("dataset", "outcome", "term", "estimate", "std.error", 
                         "statistic", "p.value", "conf.low", "conf.high")]
  data.tab$signif <- ifelse(data.tab$p.value <= alpha, 1, 0)
  
  ############## PLOTTING ###################
  
  if(inc.plot){
    
    plot.tab <- data.tab
    plot.tab$model <- factor(plot.tab$outcome, levels = dep.vars)
    keep.plot.rows <- Reduce(intersect, list(grep(paste(keep.vars.plot, collapse = "|"), plot.tab$term), 
                                grep(paste(keep.dep.plot, collapse = "|"), plot.tab$outcome), 
                                grep(paste(keep.data.plot, collapse = "|"), plot.tab$dataset)))
    plot.tab <- plot.tab[keep.plot.rows, ]
    levels(plot.tab$model) <- dep.names
    single.plot <- sum(length(dep.vars), length(data)) == 2 
    
    ### IF SINGLE DV AND SINGLE DATASET
    
    if(single.plot){
      tab.use <- plot.tab
      p <- dwplot(tab.use, intercept = F, 
                  dot_args = list(size = 5, pch = 20, 
                                  aes(colour = abs(statistic))), 
                  whisker_args = list(lwd = 1.5, 
                                      aes(colour = abs(statistic)))) + 
        scale_colour_gradient(low = "lightgray", high = "black")
      
      p <- p +  theme_bw(base_size = 12) + 
        geom_vline(xintercept = 0, linetype = "dashed") + 
        theme(axis.text.x = element_text(size = x.size), 
              axis.text.y = element_text(size = y.size),
              strip.text = element_text(size = strip.size), 
              legend.position = "none", 
              plot.margin = margin(0.75, 0.75, 0.75, 0.75, "cm")) + 
        ggtitle(paste("DV:", dep.vars))
      
      #################### LABEL CHANGES AND OTHER ADDITIONS ##################
      
      all.stand <- all(dep.vars %in% stand.vars)
      all.orig <- all(!dep.vars %in% stand.vars)
      some.mixed <- all.stand == F & all.orig == F
      
      if(all.stand){
        p <- p + xlab("\nChange in Dependent Variable (Standard Deviation Units)")}
      if(all.orig){
        p <- p + xlab("\nChange in Dependent Variable (Original Units)")}
      if(some.mixed){
        p <- p + xlab("\nChange in Dependent Variable (Mixed Units)")
        warning("Only some dependent variables are standardized: changing x-axis label to `Mixed Units`")}
      plot.list <- p
      
      if(print.plot){print(p)}# This will make sure plot shows automatically
      
    }
    
    ### IF MULTIPLE DVs OR MULTIPLE DATASETS (OR BOTH)
    
    if(!single.plot){
      
      ## If you want to produce a separate plot for each outcome
      
      if(subset.by == "outcome"){
        
        plot.list <- list()
        loop.vals <- unique(plot.tab$outcome)
        
        for(i in 1:length(loop.vals)){
          
          tab.use <- plot.tab[plot.tab$outcome == loop.vals[i], ]
          p <- dwplot(tab.use, intercept = F, 
                      dot_args = list(size = 5, pch = 20, 
                                      aes(colour = abs(statistic))), 
                      whisker_args = list(lwd = 1.5, 
                                          aes(colour = abs(statistic)))) + 
            scale_colour_gradient(low = "lightgray", high = "black") + 
            
          
          if(length(unique(tab.use$dataset)) > 1){
            p <- p + facet_wrap(~dataset, labeller = label_parsed)
          }
          
          p <- p +  theme_bw(base_size = 12) + 
            geom_vline(xintercept = 0, linetype = "dashed") + 
            theme(axis.text.x = element_text(size = x.size), 
                  axis.text.y = element_text(size = y.size),
                  strip.text = element_text(size = strip.size), 
                  legend.position = "none", 
                  plot.margin = margin(0.75, 0.75, 0.75, 0.75, "cm")) + 
            ggtitle(paste("Estimated Effects in ", loop.vals[i], "sample"))
          
          #################### LABEL CHANGES AND OTHER ADDITIONS ##################
          
          all.stand <- all(dep.vars %in% stand.vars)
          all.orig <- all(!dep.vars %in% stand.vars)
          some.mixed <- all.stand == F & all.orig == F
          
          if(all.stand){
            p <- p + xlab("\nChange in Dependent Variable (Standard Deviation Units)")}
          if(all.orig){
            p <- p + xlab("\nChange in Dependent Variable (Original Units)")}
          if(some.mixed){
            p <- p + xlab("\nChange in Dependent Variable (Mixed Units)")
            warning("Only some dependent variables are standardized: changing x-axis label to `Mixed Units`")}
          
          if(print.plot){print(p)}# This will make sure plot shows automatically
          
          plot.list[[i]] <- p
          
        } # close i loop
        
        names(plot.list) <- dep.names 
        
      } # close the multiple outcomes loop
      
      
      ## IF YOU WANT TO PRODUCE A SEPARATE PLOT FOR EACH DATASET
      
      if(subset.by == "dataset"){
        
        loop.vals <- unique(plot.tab$dataset)
        plot.list <- list()
        
        for(i in 1:length(loop.vals)){
          
          tab.use <- plot.tab[plot.tab$dataset == loop.vals[i], ]
          p <- dwplot(tab.use, intercept = F, 
                      dot_args = list(size = 5, pch = 20, 
                                      aes(colour = abs(statistic))), 
                      whisker_args = list(lwd = 1.5, 
                                          aes(colour = abs(statistic)))) + 
            scale_colour_gradient(low = "lightgray", high = "black")
          
          if(length(unique(tab.use$model)) > 1){ # Note that have to use model here and not "outcome" because model is what is taken by dwplot
            p <- p + facet_wrap(~model, labeller = label_parsed)
          }
          
          p <- p +  theme_bw(base_size = 12) + 
            geom_vline(xintercept = 0, linetype = "dashed") + 
            theme(axis.text.x = element_text(size = x.size), 
                  axis.text.y = element_text(size = y.size),
                  strip.text = element_text(size = strip.size), 
                  legend.position = "none", 
                  strip.background = element_rect(size = 0.5), 
                  plot.margin = margin(0.75, 0.75, 0.75, 0.75, "cm")) + 
            ggtitle(paste("Results for", loop.vals[i]))
          
          #################### LABEL CHANGES AND OTHER ADDITIONS ##################
          
          all.stand <- all(dep.vars %in% stand.vars)
          all.orig <- all(!dep.vars %in% stand.vars)
          some.mixed <- all.stand == F & all.orig == F
          
          if(all.stand){
            p <- p + xlab("\nEstimated Effect (DV in Standard Deviation Units)")}
          if(all.orig){
            p <- p + xlab("\nEsimated Effect (DV in Original Units)")}
          if(some.mixed){
            p <- p + xlab("\nEstimated Effect (DV in Mixed Units)")
            warning("Only some dependent variables are standardized: changing x-axis label to `Mixed Units`")}
          
          if(print.plot){print(p)}# This will make sure plot shows automatically
          
          plot.list[[i]] <- p
          
        } # close i loop
        
        names(plot.list) <- data.names 
        
      } # Close the data subset loop
    } # Close the multiple plot loop
    
    
  } # Close plot include loop
  
  
  ############# SUBSET ACCORDING TO KEEP ARGUMENTS AND BUILD FINAL OBJECT
  
  keep.tab.rows <- Reduce(intersect, list(grep(paste(keep.vars.tab, collapse = "|"), data.tab$term), 
                              grep(paste(keep.dep.tab, collapse = "|"), data.tab$outcome), 
                              grep(paste(keep.data.tab, collapse = "|"), data.tab$dataset)))
  data.tab <- data.tab[keep.tab.rows, ]
  data.tab$outcome <- factor(data.tab$outcome, levels = dep.vars)
  levels(data.tab$outcome) <- dep.names
  
  if(!is.null(stand.vars)){print(paste("The following variables have been standardized:", 
                                        paste(change.vars, collapse = ", ")))}
  if(exists("plot.list")){print(paste("NOTE: Coefficient scaling on plots is", stand.mult))}
  if(!exists("plot.list")){plot.list <- "Set inc.plot to TRUE to include plot objects in output"}
  if(!inc.mod){mod.data <- "Set inc.mod to TRUE to include model objects in output"}
  if(!all(dep.vars == dep.names)){warning("Custom dependent variable names used, check that order in dep.names matches order in dep.vars argument")}
  return(list(table = as_tibble(data.tab),
                           info = info.tab,
                           models = mod.data,
                           plot = plot.list))
  
}

